package dipl.algorithm.impl;

import dipl.algorithm.statistics.Statistics;
import dipl.algorithm.math.curve.BezierCurve;
import dipl.algorithm.math.geometry.ConvexHull2d;
import dipl.algorithm.math.geometry.op.BoundingWedge2dOps;
import dipl.algorithm.math.geometry.op.ConvexPolygons2d;
import dipl.algorithm.math.primitive.Interval1d;
import dipl.algorithm.math.primitive.Point2d;
import dipl.algorithm.math.utility.FPMath;
import dipl.algorithm.math.primitive.Tuple2;
import dipl.algorithm.exception.ArgumentNullException;
import dipl.algorithm.math.geometry.BoundingWedge2d;
import dipl.algorithm.math.geometry.op.Point2dSetOps;
import dipl.algorithm.yap.YapUtilities;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.Stack;
import java.util.Vector;
import org.apfloat.Apfloat;

/**
 * The spline intersection algorithm (with or without supercomposition)
 */
public class ScmpAlgorithmBezierTest {

	//
	// TYPES
  //

	/**
	 * Enumeration of curve segment types
	 */
	protected enum SegmentType {
		/** [t-h,t+h]-Segment */
		Left,
		/** [t,t+2h]-Segment */
		Right,
		/** base [t,t+h]-Segment */
		Base,
	}

	/**
	 * Class represents one curve segment in the subdivision tree.
	 * Along with the curve segment its convex hull and other data
	 * required to determine whether curve segments intersect in one unique point is stored.
	 */
	public class CurveSegment {

		public CurveSegment() {
			InitializeMembers();
		}

		public CurveSegment( BezierCurve curve ) {
			InitializeMembers();
			this.curve = curve;
			if( this.curve != null ) this.domain = this.curve.Domain();
		}

		public CurveSegment( CurveSegment unionSegL, CurveSegment unionSegR ) {
			InitializeMembers();
			this.isUnion = true;
			this.unionSegL = unionSegL;
			this.unionSegR = unionSegR;
		}

		public CurveSegment( BezierCurve curve, ConvexHull2d hull ) {
			InitializeMembers();
			this.curve = curve;
			if( this.curve != null ) this.domain = this.curve.Domain();
			this.hull = hull;
		}

		public CurveSegment( BezierCurve curve, ConvexHull2d hull, BoundingWedge2d wedge ) {
			InitializeMembers();
			this.curve = curve;
			if( this.curve != null ) this.domain = this.curve.Domain();
			this.hull = hull;
			this.wedge = wedge;
		}

		/** @return count of segments already checked against this segment */
		public int CountCheckedSegments() {
			//return checkedIntervals.count();
			return checkedIntervals.size();
		}

		/** @return whether or not this segment was constructed as the union of two other segments */
		public boolean IsUnion() {
			return isUnion;
		}

		/** @return left segment of segments from which this was constructed */
		public CurveSegment UnionSegL() {
			return unionSegL;
		}

		/** @return right segment of segments from which this was constructed */
		public CurveSegment UnionSegR() {
			return unionSegR;
		}

		private void InitializeMembers() {
			this.curve = null;
			this.hull = null;
			this.wedge = null;
			this.domain = new Interval1d();
			//this.comp = new IntervalComparator();
			this.checkedIntervals = new Vector<Interval1d>( 9, 9 );//new RedBlackTree<Interval1d, IntervalComparator>( comp );
			this.isUnion = false;
			this.unionSegL = null;
			this.unionSegR = null;
		}

		/**
		 * Adds given segment interval to list of intervals
		 * which have already been checked against this one by the algorithm.
		 * We store only intervals here because in the algorithm the curve segments
		 * are not always available (and must be constructed, which is expensive)
		 * @param interval
		 */
		public void AddChecked( Interval1d interval ) {
			//checkedIntervals.insertKey( o );
			checkedIntervals.add( interval );
		}

		/**
		 * @param o
		 * @return whether curve of given segment "o" is identical to this segment's curve.
		 */
		public boolean CurvesIdentical( CurveSegment o ) {
			if( this.curve == null || o.curve == null )
				return false;
			return this.curve.Equals( o.curve );
		}

		/**
		 * @param interval
		 * @return whether segment pair formed by this segment and given one has already
		 * been considered by the algorithm.
		 */
		public boolean IsAlreadyChecked( Interval1d interval ) {
			//return (checkedIntervals.search( interval ) != null);
			int len = checkedIntervals.size();
			for( int i = 0; i < len; i++ ) {
				if( checkedIntervals.contains( interval ) )
					return true;
			}
			return false;
		}

		/** @return  whether this curve segment is degenerated to one single point. */
		public boolean IsOnePoint() {
			if( curve == null )
				return false;
			return curve.IsDegenerated();
		}

		/** domain of curve segment */
		public Interval1d domain;
		/** curve segment */
		public BezierCurve curve;
		/** convex hull of curve segment */
		public ConvexHull2d hull;
		/** bounding wedge of curve segment */
		public BoundingWedge2d wedge;
		/** tree containing domains of curve segments which where already checked against this segment */
		//private RedBlackTree<Interval1d, IntervalComparator> checkedIntervals;
		private Vector<Interval1d> checkedIntervals;
		/** interval comparator */
		//private IntervalComparator comp;
		/** whether or not this segment was constructed by two others */
		private boolean isUnion;
		/** one of the segmetns from which this segment was constructed */
		private CurveSegment unionSegL;
		/** one of the segments from which this segment was constructed */
		private CurveSegment unionSegR;
	};

	/** Node of Subdivision tree for one curve. */
	protected class SubdivisionNode {

		public SubdivisionNode( CurveSegment seg ) {
			InitializeMembers();
			this.seg = seg;
		}

		public SubdivisionNode( CurveSegment seg, SubdivisionNode parent, SubdivisionNode left, SubdivisionNode right ) {
			InitializeMembers();
			this.seg = seg;
			this.parent = parent;
			this.left = left;
			this.right = right;
		}

		private void InitializeMembers() {
			this.seg = null;
			this.segSP1 = null;
			this.segSP2 = null;
			this.parent = null;
			this.left = null;
			this.right = null;
			this.depth = 0;
		}

		/** curve seg of this node */
		public CurveSegment seg;
		/** superposition curve seg [t-h,t+h] */
		public CurveSegment segSP1;
		/** superposition curve seg [t,t+2h] */
		public CurveSegment segSP2;
		/** parent node */
		public SubdivisionNode parent;
		/** left sub-curve's node */
		public SubdivisionNode left;
		/** right sub-curve's node */
		public SubdivisionNode right;
		/** Depth of node in subdivision tree. Depth 0 means root */
		public int depth;
	};

	//
	// CONSTRUCTOR METHODS
	//

	public ScmpAlgorithmBezierTest() {
		statistics = null;
		fpMath = null;
		intersections = null;
		macroPairs = null;
		microPairs = null;
	}

	//
	// PROPERTIES
	//

	/**
	 * @return intersections found
	 */
	public Iterator<Intersection> Intersections() {
		return intersections.iterator();
	}

	/**
	 * @return object collecting statistic data during algorithm
	 */
	public Statistics Statistics() {
		return statistics;
	}

	//
	// PUBLIC METHODS
	//

	/**
	 * Returns all intersections between given two curves f and g.
	 * It is assumed that both curves are relatively prime
	 * (have no common segment and contain no translated segments of each other)!
	 *
	 * @param f first input curve
	 * @param g second input curve
	 * @param parameters flags (from {@link Parameters}) controlling algorithm behaviour
	 * @param treshold the algorithm stops if curve segment "size" run below this value
	 * @param L L-constant. If Yap-method is used it indicates that coefficients of input curves are (L,L)-bit floats.
	 * @throws ArgumentNullException one of given curves is <code>null</code>
	 */
	public void CalculateIntersections(
		BezierCurve f,
		BezierCurve g,
		int parameters,
		Apfloat treshold,
		int L
	) throws ArgumentNullException {

		if( f == null || g == null )
			throw new ArgumentNullException();

		this.parameters = parameters;
		this.startPrecision = 48; // 16 decimal places in dual number system
		this.currentPrecision = startPrecision;
		if( (parameters & Parameters.COLLECT_DATA) != 0 ) {
			// start collection of data
			this.statistics = new Statistics();
			this.statistics.startPrecision = currentPrecision;
			this.statistics.startTimeMS = System.currentTimeMillis();
		}

		this.macroPairs = new LinkedList<Tuple2<Interval1d, Interval1d>>();
		if( (parameters & Parameters.BOUNDS_YAP) != 0 ) {
			this.microPairs = new LinkedList<Tuple2<Interval1d, Interval1d>>();
		}
		this.intersections = new LinkedList<Intersection>();

		this.f = f;
		this.g = g;

		// check whether curves are degenerated
		if( f.IsDegenerated() && g.IsDegenerated() ) {
			if( f.Controls()[0].Equals( g.Controls()[0] ) ) {
				AddIntersection( f.Domain(), g.Domain(), IntersectionType.Unique );
			}
			if( (parameters & Parameters.COLLECT_DATA) != 0 ) {
				this.statistics.endTimeMS = System.currentTimeMillis();
			}
			return;
		}

		this.maxDepth_f = 0;
		this.maxDepth_g = 0;
		this.fpMath = new FPMath( startPrecision );

		// Current segment pair (current candidate)
		Tuple2<SubdivisionNode,SubdivisionNode> cand = null;
		// stack containing candidate segment pairs
		Stack<Tuple2<SubdivisionNode,SubdivisionNode>> candidates = new Stack<Tuple2<SubdivisionNode,SubdivisionNode>>();

		// create first candidate pair
		SubdivisionNode fRoot = new SubdivisionNode( new CurveSegment( f ) );
		SubdivisionNode gRoot = new SubdivisionNode( new CurveSegment( g ) );
		candidates.push( new Tuple2<SubdivisionNode,SubdivisionNode>( fRoot, gRoot ) );

		Apfloat squareTreshold =  treshold.multiply( treshold );
		Apfloat size_f = null; // "size" of segment of f
		Apfloat size_g = null; // "size" of segment of g
		Apfloat delta = null;
		// calculate yap separation bounds ?
		if( (parameters & Parameters.BOUNDS_YAP) != 0 ) {
			delta = YapUtilities.CalculateMaxSquareDelta( L, f.Degree(), g.Degree() );
		}

		while( candidates.size() > 0 ) {
			cand = candidates.pop();

			// "size" of segments
			ConstructConvexHull( cand.fst.seg );
			ConstructConvexHull( cand.snd.seg );
			size_f = ConvexPolygons2d.CalculateDiameter( fpMath, cand.fst.seg.hull.Vertices(), true );
			size_g = ConvexPolygons2d.CalculateDiameter( fpMath, cand.snd.seg.hull.Vertices(), true );

			// possible intersection because at least one segment is too short?
			if( (parameters & Parameters.BOUNDS_YAP) == 0 ) {
				if( size_f.compareTo( squareTreshold ) < 0 && size_g.compareTo( squareTreshold ) < 0 ) {
					AddIntersection( cand.fst.seg.domain, cand.snd.seg.domain, IntersectionType.Possible );
					continue;
				}
			}
			else {
				if( size_f.compareTo( delta ) < 0 && size_g.compareTo( delta ) < 0 ) {
					AddMicroPair( cand.fst.seg.domain, cand.snd.seg.domain );
				}
			}

			// supercomposition?
			if( (parameters & Parameters.SUPERCOMPOSITION) != 0 ) {
				// before subidividing pair we test supercompositioned segment pairs.

				boolean foundInters = false; // found intersection?
				boolean foundDisjoint = false; // whether or not one pair was disjoint
				int i = 0;
				for( i = 0; i < 9; i++ ) {
					CurveSegment scmpFst = CalculateSCMPSegment( cand.fst, f, scmpCalculationSequence[i][0] );
					if( scmpFst == null ) // first sp-segment is null -> nothing to do
						continue;
					CurveSegment scmpSnd = CalculateSCMPSegment( cand.snd, g, scmpCalculationSequence[i][1] );
					// second sp-segment is null or segments already checked -> nothing to do
					if( scmpSnd == null || IsSegmentPairAlreadyChecked( scmpFst, scmpSnd ) )
						continue;

					scmpFst.AddChecked( scmpSnd.domain );
					scmpSnd.AddChecked( scmpFst.domain );

					int test = TestPair( scmpFst, scmpSnd );
					// we test given pair for no intersection
					if( test == 0 ) { // condition (A)
            foundDisjoint = true;
						break;
					}
					// then for unique intersection
					if( test == 1 ) { // condition (B)
						AddMacroPair( scmpFst.domain, scmpSnd.domain );
						foundInters = true;
						break;
					}
				}
				// if a unique intersection was found or a pair was disjoint, no further subdivision is needed
				if( foundInters || foundDisjoint )
					continue;
			}
			// no supercomposition -> only check candidate pair
			else {
				int test = TestPair( cand.fst.seg, cand.snd.seg );
				if( test == 1 ) {
					AddMacroPair( cand.fst.seg.domain, cand.snd.seg.domain );
					continue;
				}
				else if( test == 0 ) {
					continue;
				}
			}

			// no decision was made
			// -> subdivide  curve with larger size
			if( size_f.compareTo( size_g ) > 0 ) {
				// f already subdivided?
				if( cand.fst.left == null ) {
					SubdivideNode( cand.fst, true );
				}
				// add new canditates (f1,g), (f2,g)
				candidates.push( new Tuple2<SubdivisionNode,SubdivisionNode>( cand.fst.left, cand.snd ) );
				candidates.push( new Tuple2<SubdivisionNode,SubdivisionNode>( cand.fst.right, cand.snd ) );
			}
			else {
				// g already partitioned?
				if( cand.snd.left == null ) {
					SubdivideNode( cand.snd, false );
				}
				// add new canditates (f,g1), (f,g2)
				candidates.push( new Tuple2<SubdivisionNode,SubdivisionNode>( cand.fst, cand.snd.left ) );
				candidates.push( new Tuple2<SubdivisionNode,SubdivisionNode>( cand.fst, cand.snd.right ) );
			}
		}
		// collect intersections found
		for( Tuple2<Interval1d,Interval1d> p : macroPairs ) {
			AddIntersection( p.fst, p.snd, IntersectionType.Unique );
		}
		if( (parameters & Parameters.BOUNDS_YAP) != 0 ) {
			for( Tuple2<Interval1d,Interval1d> p : microPairs ) {
				AddIntersection( p.fst, p.snd, IntersectionType.Micro );
			}
		}
		// statistics?
		if( (parameters & Parameters.COLLECT_DATA) != 0 ) {
			statistics.endTimeMS = System.currentTimeMillis();
			statistics.subdivisionDepth_f = maxDepth_f;
			statistics.subdivisionDepth_g = maxDepth_g;
			statistics.endPrecision = currentPrecision;
			CalculateSubdivsionTreeMeans( fRoot, gRoot );
		}
	}

	//
	// AUXILLIARY METHODS
	//

	/**
	 * Adds given intersection to intersection list.
	 * @param f subinterval of first input curve
	 * @param g subinterval of second input curve
	 * @param type
	 */
	protected void AddIntersection( Interval1d f, Interval1d g, IntersectionType type ) {
		this.intersections.addLast(
			new Intersection( new Interval1d( f ), new Interval1d( g ), type ) );
	}

	/**
	 * Adds pair to macro pair list,
	 * if a pair containing this pair's intervals is not already contained.
	 * @param f first segment interval
	 * @param g second segment interval
	 */
	protected void AddMacroPair( Interval1d f, Interval1d g ) {
		if( (parameters & Parameters.SUPERCOMPOSITION) != 0 && macroPairs.size() > 0 ) {
			for( Tuple2<Interval1d,Interval1d> p : this.macroPairs ) {
				if( p.fst.ProperlyIntersects( f ) &&
						p.snd.ProperlyIntersects( g ) )
					return;
			}
		}
		this.macroPairs.addLast( new Tuple2<Interval1d,Interval1d>( f, g ) );
	}

	/**
	 * Adds pair to micro pair list,
	 * if a pair containing this pair's intervals is not already contained.
	 * @param f first segment interval
	 * @param g second segment interval
	 */
	protected void AddMicroPair( Interval1d f, Interval1d g ) {
		if( (parameters & Parameters.SUPERCOMPOSITION) != 0 && microPairs.size() > 0 ) {
			for( Tuple2<Interval1d,Interval1d> p : this.microPairs ) {
				if( p.fst.ProperlyIntersects( f ) &&
						p.snd.ProperlyIntersects( g ) )
					return;
			}
		}
		this.microPairs.addLast( new Tuple2<Interval1d,Interval1d>( f, g ) );
	}

	/**
	 * Trys to add supercomposition pair formed by given two curve segments to given collection.
	 * If both segments exist and segment pair was not already checked the pair is added to given collection.
	 * @param f
	 * @param g
	 * @param pairs
	 */
	protected void AddSupercompositionPair( CurveSegment f, CurveSegment g, LinkedList<Tuple2<CurveSegment,CurveSegment>> pairs ) {
		if( f != null && g != null && !IsSegmentPairAlreadyChecked( f, g ) ) {
			pairs.addFirst( new Tuple2<CurveSegment,CurveSegment>( f, g ) );
			f.AddChecked( g.domain );
			g.AddChecked( f.domain );
		}
	}

  /**
	 * Calculates left ([t-h,t+h]) superposition segment of given subdivision node (using extrapolation) if it
	 * doesn't already exist.
	 * @param n
	 * @param initialDomain
	 */
	protected void CalculateLeftSPSegmentExtrapolated( SubdivisionNode n, Interval1d initialDomain ) {
		if( n.parent != null ) { // if parent is null then we are at the root
			if( n.segSP1 == null ) { // -> parent is first sp-segment
				if( n.parent.right == n ) { // we are right child
					n.segSP1 = n.parent.seg;
				}
				else { // we are left child
					Interval1d dom = new Interval1d( n.seg.domain ); // domain of curve segment
					if( initialDomain.ContainsValue( fpMath.m.subtract( dom.L(), dom.Length() ) ) ) {
						dom.SetBoundsUnsafe( fpMath.m.subtract( dom.L(), dom.Length() ), dom.U() );
						n.segSP1 = new CurveSegment( n.seg.curve.ChangeDomain( fpMath, dom ) );
					}
				}
			}
		}
	}

	/**
	 * Calculates right ([t,t+2h]) superposition segment of given subdivision node (using extrapolation) if it
	 * doesn't already exist.
	 * don't already exist.
	 * @param n
	 * @param initialDomain
	 */
	protected void CalculateRightSPSegmentExtrapolated( SubdivisionNode n, Interval1d initialDomain ) {
		if( n.parent != null ) { // if parent is null then we are at the root
			if( n.segSP2 == null ) { // -> extrapolate second sp-segment if its contained in initial domain
				if( n.parent.right == n ) { // we are right child
					Interval1d dom = new Interval1d( n.seg.domain ); // domain of curve segment
					if( initialDomain.ContainsValue( fpMath.m.add( dom.U(), dom.Length() ) ) ) {
						dom.SetBoundsUnsafe( dom.L(), fpMath.m.add( dom.U(), dom.Length() ) );
						n.segSP2 = new CurveSegment( n.seg.curve.ChangeDomain( fpMath, dom ) );
					}
				}
				else { // we are left child
					n.segSP2 = n.parent.seg;
				}
			}
		}
	}

	/**
	 * Calculates left ([t-h,t+h]) superposition curve segment of given subdivision node
	 * (without using extrapolation) if it doesn't already exist
	 * @param n
	 * @param originalDomain domain of orignal curve
	 * @param firstCurve whether or not node n belongs to sd-tree of first input curve
	 */
	protected void CalculateLeftSPSegmentNotExtrapolated( SubdivisionNode n, Interval1d originalDomain, boolean firstCurve ) {
		if( n.parent != null ) {
			// now aquire and construct superposition data
			if( n.segSP1 == null ) {
				if( n.parent.right == n ) {
					n.segSP1 = n.parent.seg;
				}
				else {
					SubdivisionNode sibling = GetOrCreateLeftSibling( n, originalDomain, firstCurve );
					if( sibling != null ) {
						n.segSP1 = new CurveSegment( sibling.seg, n.seg );
						n.segSP1.domain = new Interval1d( sibling.seg.domain.L(), n.seg.domain.U() );
						// we store superpositioned hull in  sibling found so that
						// it don't have to be constructed when sibling is to be superpositioned
						sibling.segSP2 = n.segSP1;
					}
				}
			}
		}
	}

	/**
	 * Calculates right ([t,t+2h]) superposition curve segment of given subdivision node
	 * (without using extrapolation) if it doesn't already exist
	 * @param n
	 * @param originalDomain domain of orignal curve
	 * @param firstCurve whether or not node n belongs to sd-tree of first input curve
	 */
	protected void CalculateRightSPSegmentNotExtrapolated( SubdivisionNode n, Interval1d originalDomain, boolean firstCurve ) {
		if( n.parent != null ) {
			if( n.segSP2 == null ) {
				// now aquire and construct superposition data
				if( n.parent.right == n ) {
					SubdivisionNode sibling = GetOrCreateRightSibling( n, originalDomain, firstCurve );
					if( sibling != null ) {
						n.segSP2 = new CurveSegment( n.seg, sibling.seg );
						n.segSP2.domain = new Interval1d( n.seg.domain.L(), sibling.seg.domain.U() );
						// we store superpositioned segment in sibling found so that
						// it doesn't have to be constructed when sibling is to be superpositioned
						sibling.segSP1 = n.segSP2;
					}
				}
				else {
					n.segSP2 = n.parent.seg;
				}
			}
		}
	}

	/**
	 * Calculates and returns supercomposition curve segment of given type for segment in given subdivision node.
	 * @param node
	 * @param origCurve
	 * @param type
	 */
	protected CurveSegment CalculateSCMPSegment( SubdivisionNode node, BezierCurve origCurve, SegmentType type ) {
		switch( type ) {
			case Left:
				if( node.seg.curve.CanBeExtrapolated() ) {
					CalculateLeftSPSegmentExtrapolated( node, origCurve.Domain() );
				}
				else {
					CalculateLeftSPSegmentNotExtrapolated( node, origCurve.Domain(), origCurve == this.f );
				}
				return node.segSP1;
			case Right:
				if( node.seg.curve.CanBeExtrapolated() ) {
					CalculateRightSPSegmentExtrapolated( node, origCurve.Domain() );
				}
				else {
					CalculateRightSPSegmentNotExtrapolated( node, origCurve.Domain(), origCurve == this.f );
				}
				return node.segSP2;
		}
		return node.seg;
	}

	/**
	 * Returns the convex hull from the union of control polygons of given two sibling bezier curve
	 * segments c1 and c2.
	 * It is important to note that resulting polygon will automatically contain the indices of first
	 * control point of c1 and last control point of c2 on the hull if they are contained in hull's vertices.
	 * This is possible because of end-point properties of bezier-curves.
	 * Those indices are most important to determine whether convex hulls of curves are "crossing" later on.
	 *
	 * Method assumes hulls of given segments to be not null!
	 * Running time: O(n+m)
	 * @param c1
	 * @param c2
	 * @return
	 */
	protected ConvexHull2d CombineSegmentHulls( CurveSegment s1, CurveSegment s2 ) {
		ConvexHull2d res = new ConvexHull2d( fpMath, Point2dSetOps.Concat2( s1.hull.Vertices(), s2.hull.Vertices() ) );
		// find end-point vertices
		res.SetFirstPointsetIndex( -1 );
		res.SetLastPointsetIndex( -1 );
		Point2d[] resVert = res.Vertices();
		Point2d p10 = s1.hull.Pointset()[0];
		Point2d p2n = s2.hull.Pointset()[s2.hull.Pointset().length-1];
		for( int i = 0; i < resVert.length; i++ ) {
			if( resVert[i] == p10 ) {
				res.SetFirstPointsetIndex( i );
			}
			if( resVert[i] == p2n ) {
				res.SetLastPointsetIndex( i );
			}
		}
		return res;
	}

	protected void ConstructConvexHull( CurveSegment seg ) {
		if( seg.hull != null ) {
			return;
		}
		if( seg.IsUnion() ) {
			// combine hulls
			ConstructConvexHull( seg.UnionSegL() );
			ConstructConvexHull( seg.UnionSegR() );
			seg.hull = CombineSegmentHulls( seg.UnionSegL(), seg.UnionSegR() );
		}
		else {
			seg.hull = new ConvexHull2d( fpMath, seg.curve.Controls() );
		}
	}

	/**
	 * Tries to construct bounding wedge for given curve segment if it was not
	 * already constructed.
	 *
	 * @param seg
	 */
	protected void ConstructBoundingWedge( CurveSegment seg ) {
		if( seg.wedge != null ) {
			return;
		}
		// wedge is union of segments? (applies only if curve-segment was not extrapolated)
		if( seg.IsUnion() ) {
			// -> now combine these wedges (construct from union of hodograph vectors)
			//
			// INSERT SOME CURVE-SPECIFIC CODE HERE
			//
			// this code here is only for testing purposes since Bezier-Curves can be extrapolated
			// and this control-path is never executed!
			Point2d[] controls = Point2dSetOps.Concat2( seg.UnionSegL().curve.Controls(), seg.UnionSegR().curve.Controls() );
			seg.wedge = BezierCurve.BoundingWedge( fpMath, controls );
		}
		else {
			seg.wedge = seg.curve.BoundingWedge( fpMath );
		}
	}

	/**
	 * Returns left sibling node of given one. If domain of sibling will
	 * exceed given domain, <code>null</code> is returned.
	 * If no left sibling node exists in subdivision tree, it is created together
	 * with its parent node. (The grandparent has already be created)
	 * WARNING:
	 * If given node is not the root it is assumed that it is the left child of its parent.
	 * @param node
	 * @param domain
	 * @param firstCurve whether or not the search takes place in subd.-tree of first curve
	 * @return
	 */
	protected SubdivisionNode GetOrCreateLeftSibling( SubdivisionNode node, Interval1d domain, boolean firstCurve ) {
		Interval1d dom = node.seg.domain; // domain of curve segment
		// check whether domain of right sibling is contained in
		// intial curve domain.
		if( !domain.ContainsValue( fpMath.m.subtract( dom.L(), dom.Length() ) ) )
			return null;
		// now that we know that left sibling must be contained in tree, it is
		// contained in a subtree left of given node
		SubdivisionNode p = node.parent; // current node under consideration
		int depth = 1; // depth of given node (starts with 1 here because given node cannot be root)
		while( p.parent.left == p ) {
			p = p.parent;
			depth++;
		}
		depth++; // make depth one based
		// go right down the left sub-tree and eventually create a subtree on the way.
		// A left sub-tree must exist here for two reasons:
		// Firstly if not the left sibling would be already outside the initial domain.
		// Secondly it must have been created by a previous subdivision of p.
		p = p.parent.left;
		while( depth > 1 ) {
			if( p.right == null ) { // create subtree?
				SubdivideNode( p, firstCurve );
			}
			p = p.right;
			depth--;
		}
		return p;
	}

	/**
	 * Returns right sibling node of given one. If domain of sibling would
	 * exceed given domain, <code>null</code> is returned.
	 * If no right sibling node exists in subdivision tree, it is created together
	 * with its parent node. (The grandparent has already be created)
	 * WARNING:
	 * If given node is not the root it is assumed that it is the right child of its parent.
	 * @param node
	 * @param domain
	 * @param firstCurve whether or not the search takes place in subd.-tree of first curve
	 * @return
	 */
	protected SubdivisionNode GetOrCreateRightSibling( SubdivisionNode node, Interval1d domain, boolean firstCurve ) {
		Interval1d dom = node.seg.domain; // domain of curve segment
		// check whether domain of right sibling is contained in
		// intial curve domain.
		if( !domain.ContainsValue( fpMath.m.add( dom.U(), dom.Length() ) ) )
			return null;
		// now that we know that right sibling must be contained in tree, it is
		// contained in a subtree right of given node
		SubdivisionNode p = node.parent; // current node under consideration
		int depth = 1; // depth of given node (starts with 1 here because given node cannot be root)
		while( p.parent.right == p ) {
			p = p.parent;
			depth++;
		}
		depth++; // make depth one-based
		// go left down the right sub-tree and eventually create subtree on the way
		// a right sub-tree must exist here otherwise the right sibling would
		// be already outside the initial domain.
		p = p.parent.right; // right sub-tree must exist here because it was create by subdivision of p
		while( depth > 1 ) {
			if( p.left == null ) { // create subtree?
				SubdivideNode( p, firstCurve );
			}
			p = p.left;
			depth--;
		}
		return p;
	}

	protected boolean IsSegmentPairAlreadyChecked( CurveSegment s1, CurveSegment s2 ) {
		return (s1.IsAlreadyChecked( s2.domain ) || s2.IsAlreadyChecked( s1.domain ));
	}

	/**
	 * Returns whether given two curve segments are guaranteed to have no intersection
	 * @param a
	 * @param b
	 * @return
	 */
	protected boolean NoIntersection( CurveSegment a, CurveSegment b ) {
		ConstructConvexHull( a );
		ConstructConvexHull( b );
		return !ConvexPolygons2d.Intersect( fpMath, a.hull.Vertices(), b.hull.Vertices() );
	}

	/**
	 * Subdivides given node of subdivision tree.
	 * @param node
	 * @param firstCurve whether or not node belongs to subdivision tree of first curve
	 */
	protected void SubdivideNode( SubdivisionNode node, boolean firstCurve ) {
		// keep track of maximal subdivision depth
		// eventually increase precision by 1 (one bit per subdivision) if a new
		// subdivision level was reached
		if( firstCurve ) {
			maxDepth_f = Math.max( node.depth+1, maxDepth_f );
			if( currentPrecision-startPrecision < maxDepth_f ) {
				currentPrecision++;
				this.fpMath.SetPrecision( currentPrecision );
			}
			if( (parameters & Parameters.COLLECT_DATA) != 0 ) {
				statistics.numSubdivisions_f++;
			}
		}
		else {
			maxDepth_g = Math.max( node.depth+1, maxDepth_g );
			if( currentPrecision-startPrecision < maxDepth_g ) {
				currentPrecision++;
				this.fpMath.SetPrecision( currentPrecision );
			}
			if( (parameters & Parameters.COLLECT_DATA) != 0 ) {
				statistics.numSubdivisions_g++;
			}
		}
		Tuple2<BezierCurve,BezierCurve> subc = null;
		try {
			 subc = node.seg.curve.Subdivide( fpMath, fpMath._0_5 );
		}
		catch( Exception e ) {
			e.printStackTrace();
			return;
		}
		node.left = new SubdivisionNode( new CurveSegment( subc.fst ), node, null, null );
		node.right = new SubdivisionNode( new CurveSegment( subc.snd ), node, null, null );
		node.left.depth = node.right.depth = node.depth+1;
	}

	/**
	 * Returns whether given segment pair
	 * (A) has no intersection (0 is returned)
	 * (B) has a unique intersection (1 is returned)
	 * If neither (A) nor (B) -1 is returned.
	 * @param a
	 * @param b
	 * @return
	 */
	protected int TestPair( CurveSegment a, CurveSegment b ) {
		// we test given pair for no intersection
		if( NoIntersection( a, b ) ) { // condition (A)
			return 0;
		}
		// then for unique intersection
		if( UniqueIntersection( a, b ) ) { // condition (B)
			return 1;
		}
		return -1;
	}

	/**
	 * Returns whether given two curve segments have a unique intersection
	 * @param a
	 * @param b
	 * @return
	 */
	protected boolean UniqueIntersection( CurveSegment a, CurveSegment b ) {
		// crossing ?
		try {
			if( !ConvexPolygons2d.Crossing( a.hull, b.hull ) ) {
				return false;
			}
		}
		catch( Exception e ) {
			e.printStackTrace();
			return false;
		}
		// tangent wedges disjoint ?
		ConstructBoundingWedge( a );
		ConstructBoundingWedge( b );
		return !BoundingWedge2dOps.Overlap( fpMath, a.wedge, b.wedge );
	}


	//
	// STATICSTICAL METHODS
	//

	/**
	 * Method calculates several mean values based on subdivision trees
	 * of current input curves.
	 * The values are stored in the statistics data object
	 * @param fRoot root of subdivsion tree of first curve
	 * @param gRoot root of subdivsion tree of second curve
	 */
	protected void CalculateSubdivsionTreeMeans( SubdivisionNode fRoot, SubdivisionNode gRoot ) {
		// traverse subdivision trees in dfs order
		Stack<SubdivisionNode> st = new Stack<SubdivisionNode>();
		SubdivisionNode n = null;
		// for f
		st.push( fRoot );
		statistics.minComparedIntervals_f = Integer.MAX_VALUE;
		statistics.maxComparedIntervals_f = 0;
		statistics.numNodes_f = 0;
		double sum = 0.0;
		int numSegments = 0;
		while( st.size() > 0 ) {
			n = st.pop();
			numSegments = n.seg.CountCheckedSegments();
			statistics.numNodes_f++;
			sum += numSegments;
			statistics.minComparedIntervals_f = Math.min( statistics.minComparedIntervals_f, numSegments );
			statistics.maxComparedIntervals_f = Math.max( statistics.maxComparedIntervals_f, numSegments );
			if( n.left != null ) st.push( n.left );
			if( n.right != null ) st.push( n.right );
		}
		statistics.meanComparedIntervals_f = sum/(double)statistics.numNodes_f;
		// for g
		st.push( gRoot );
		statistics.minComparedIntervals_g = Integer.MAX_VALUE;
		statistics.maxComparedIntervals_g = 0;
		statistics.numNodes_g = 0;
		sum = 0.0;
		while( st.size() > 0 ) {
			n = st.pop();
			numSegments = n.seg.CountCheckedSegments();
			statistics.numNodes_g++;
			sum += numSegments;
			statistics.minComparedIntervals_g = Math.min( statistics.minComparedIntervals_g, numSegments );
			statistics.maxComparedIntervals_g = Math.max( statistics.maxComparedIntervals_g, numSegments );
			if( n.left != null ) st.push( n.left );
			if( n.right != null ) st.push( n.right );
		}
		statistics.meanComparedIntervals_g = sum/(double)statistics.numNodes_g;
	}

	//
	// STATIC MEMBERS
	//

	/**
	 * Sequence of Supercomposition-Segments to calculate
	 */
	private static final SegmentType[][] scmpCalculationSequence = new SegmentType[][] {
		{ SegmentType.Left,	 SegmentType.Base },	// (f[t-h,t+h],g)
		{ SegmentType.Right, SegmentType.Base },	// (f[t,t+2h], g)
		{ SegmentType.Right, SegmentType.Right },	// (f[t,t+2h], g[t,t+2h])
		{ SegmentType.Right, SegmentType.Left },	// (f[t,t+2h], g[t-h,t+h])
		{ SegmentType.Base,	 SegmentType.Left },	// (f,         g[t-h,t+h])
		{ SegmentType.Base,	 SegmentType.Right },	// (f,         g[t,t+2h])
		{ SegmentType.Left,	 SegmentType.Right },	// (f[t-h,t+h],g[t,t+2h])
		{ SegmentType.Left,	 SegmentType.Left },	// (f[t-h,t+h],g[t-h,t+h])
		{ SegmentType.Base,	 SegmentType.Base },	// original pair
	};

	//
	// MEMBERS
	//

	/**
	 * one of the input curves
	 */
	private BezierCurve f;
	/**
	 * one of the input curves
	 */
	private BezierCurve g;
	/**
	 * Object supporting fixed-precision floating point arithmetic
	 */
	private FPMath fpMath;
	/**
	 * collects data during current algorithm run
	 */
	private Statistics statistics;
	/**
	 * parameters controlling algorithm behaviour
	 */
	private int parameters;
	/**
	 * starting precision
	 */
	private int startPrecision;
	/**
	 * current precision
	 */
	private int currentPrecision;
	/**
	 * maximum depth of a subdivision node for first curve
	 */
	private int maxDepth_f;
	/**
	 * maximum depth of a subdivision node for second curve
	 */
	private int maxDepth_g;
	/**
	 * macro-pairs found
	 */
	private LinkedList<Tuple2<Interval1d,Interval1d>> macroPairs;
	/**
	 * micro-pairs found
	 */
	private LinkedList<Tuple2<Interval1d,Interval1d>> microPairs;
	/**
	 * intersections found
	 */
	private LinkedList<Intersection> intersections;
}
